﻿//
// Created by 17616 on 2021/3/19.
//
#include <vector>
#include <atomic>
#include <string>
#include <exception>
#include <iostream>
#include "public.h"
#include <opencv2/opencv.hpp>
#include <thread>
#include <Windows.h>
#include <fstream>
// [2021年1月2日164039][VP101M F4 2^1us][内反射暗场 2A电流][玻璃片0]
// [2021年1月2日164039][VP101M F4 2123us][内反射暗场 2A电流][玻璃片0]



//生成文件列表
//生成时间列表
//将所有图片载入内存
//检查所有图片的size and type
//根据size分成多块ROI，存储在无锁队列中
//创建线程池,
//每个线程对索取ROI进行处理
//生成：
// 1.tiff 灰度
// 2.欠曝-合适-过曝.png 灰度
// 3.有效张数.png
// 4.恢复时和原图的差异。
#define BLOCK_SIZE (10)
#define ExposureWeight (2.02)
std::vector<cv::Mat> mats;
std::vector<cv::Rect> ROIs(BLOCK_SIZE * BLOCK_SIZE);
std::vector<double> expTime;
cv::Mat tiff;//等效曝光时间为1us
cv::Mat expState; //保存过曝欠曝信息
cv::Mat effectiveness; //保存有效图片数量x10
std::atomic_uint64_t histogram[512]; //存储误差的直方图
void histogram_init() {
    for (int i = 0; i < 512; ++i) {
        histogram[i] = 0;
    }
}
void calculate_roi(int rows, int cols);

void calculate_tiff(int idx, const std::vector<std::string> &args);

int m_main(const std::vector<std::string> &args) {
    histogram_init();
    //获取曝光时间
    for (auto &i:args) {
        //std::cout << GetTimer(i) << std::endl;
        expTime.push_back(GetTimer(GetFileName(i)));
    }

    //载入图片到内存
    printf("//载入图片到内存...\r\n");
    mats.resize(args.size());
    thread_pool tp(0, args.size());
    tp.post([&](int start, int end) {
        for (int num = start; num < end; ++num) {
            mats[num] = cv::imread(args[num], cv::IMREAD_GRAYSCALE);
            //如果发现有错误提前结束线程
            if (mats[num].cols == 0 || mats[num].rows == 0) { return; }
            printf("Load(%fus) %s\r\n", expTime[num], GetFileName(args[num]).c_str());
        }
    });
    tp.join();//等待线程池执行完毕
    printf("//载入图片到内存 √\r\n");
    //检查所有图片的size and type
    printf("//检查所有图片的size and type...");
    if (mats[0].cols == 0 || mats[0].rows == 0) { return -1; }
    for (int i = 1; i < mats.size(); ++i) {
        if (mats[i].size != mats[0].size ||
            mats[i].type() != mats[0].type()) { return -1; }
    }
    printf("\t√\r\n");



    //准备主程序的变量
    tiff = cv::Mat(mats[0].size(), CV_32F);
    expState = cv::Mat(mats[0].size(), CV_8U);
    effectiveness = cv::Mat(mats[0].size(), CV_8U);

    //主程序，生成
    printf("//主程序，生成\r\n");
    printf("Hardware concurrency=%d\r\n", std::thread::hardware_concurrency());

    //分为不同的行用于并行处理
    //每个线程都需要的变量
    volatile size_t done_lines = 0;//处理好的行数
    std::vector<double> logTime(args.size());
    int maxExpTimeIdx = 0;
    double maxExpTime = 0;
    int minExpTimeIdx = 0;
    double minExpTime = DBL_MAX;
    for (int j = 0; j < args.size(); ++j) {
        logTime[j] = log(expTime[j]);
        if (maxExpTime <= expTime[j]) {
            maxExpTime = expTime[j];
            maxExpTimeIdx = j;
        }
        if (minExpTime > expTime[j]) {
            minExpTime = expTime[j];
            minExpTimeIdx = j;
        }
    }

    thread_pool tp2(0, mats[0].rows);
    tp2.post([&](int start, int end) {
        double weightSum;
        double ASum;
        std::vector<double> weight(args.size());
        std::vector<double> logPixel(args.size());
        std::vector<const uchar *> data(args.size());
        std::vector<int> step(args.size());
        for (int j = 0; j < args.size(); ++j) {
            data[j] = mats[j].data;
            step[j] = mats[j].step;
        }

        uchar pixel;
        for (int r = start; r < end; ++r) {
            for (int c = 0; c < mats[0].cols; ++c) {
                weightSum = 0;
                ASum = 0;
                for (int p = 0; p < args.size();//遍历每张图
                     ++p) {
                    pixel = *(data[p] + r * step[p] + c);
                    weight[p] = WTable[pixel];
                    weightSum += WTable[pixel];
                    logPixel[p] = LogTable[pixel];
                    ASum += (logPixel[p] - logTime[p]) * WTable[pixel];
                }
                size_t _effectiveness_ = weightSum * 10.0;
                effectiveness.at<uchar>(r, c) = _effectiveness_ > 255 ? 255 : _effectiveness_;
                if (weightSum > ExposureWeight) {//至少有ExposureWeight张图完全有效
//                            *(float *) (tiff.at<cv::Vec4b>(r, c).val) =  exp(ASum / weightSum);
                    tiff.at<float>(r, c) = exp(ASum / weightSum);
                    expState.at<uchar>(r, c) = 128;//曝光合适
                } else {
                    auto minPixel = *(data[minExpTimeIdx] + r * step[minExpTimeIdx] + c);
                    auto maxPixel = *(data[maxExpTimeIdx] + r * step[maxExpTimeIdx] + c);
                    if (minPixel >= 100) {
                        //曝光过度
//                                *(float *) (tiff.at<cv::Vec4b>(r, c).val) = minPixel / minExpTime;
                        tiff.at<float>(r, c) = minPixel / minExpTime;
                        expState.at<uchar>(r, c) = 255;
                    } else if (maxPixel <= 64) {
                        //欠曝
//                                *(float *) (tiff.at<cv::Vec4b>(r, c).val) = maxPixel / maxExpTime;
                        tiff.at<float>(r, c) = maxPixel / maxExpTime;
                        expState.at<uchar>(r, c) = 0;
                    } else {
                        //一部分图片欠曝，一部分图片过曝，没有足够曝光合适
                        printf("R,C = %d,%d,%.3f<%.3f\r\n", r, c, weightSum, ExposureWeight);
                        //正常处理
//                                *(float *) (tiff.at<cv::Vec4b>(r, c).val) = exp(ASum / weightSum);
                        tiff.at<float>(r, c) = exp(ASum / weightSum);
                        expState.at<uchar>(r, c) = 128 + 64;
                    }
                }
                //以上处理完毕，开始计算误差
                for (int p = 0; p < args.size();//遍历每张图
                     ++p) {
                    pixel = *(data[p] + r * step[p] + c);
                    int buf = 0.5 + expTime[p] * tiff.at<float>(r, c);
                    uchar recover_pixel = buf > 255 ? 255 : buf;
                    int res = recover_pixel - pixel;
                    if (res < 0) { res += 512; }
                    histogram[res]++;
                }

            }

            ++done_lines;//处理好一行了
        }
    });
    std::thread report_td([&]() {
        float current_rate = 0;
        float last_rate = 0;
        while (current_rate < 0.98) {
            current_rate = (float) done_lines / mats[0].rows;
            if (current_rate - last_rate > 0.02) {
                //播报一次
                printf("%.1f%%...\r\n", current_rate * 100);
                last_rate = current_rate;
            }
            Sleep(1000);
        }
    });//用于播报进度的进程
    report_td.join();
    tp2.join();


    //最后是保存相应的东西到文件。
    printf("//保存文件。\r\n");
    std::thread td1([&]() { cv::imwrite(GetFilePath(args[0]) + "合成.tiff", tiff); });
    std::thread td2([&]() { cv::imwrite(GetFilePath(args[0]) + "ExposureState.png", expState); });
    std::thread td3([&]() { cv::imwrite(GetFilePath(args[0]) + "effectiveness_x10.png", effectiveness); });
    std::thread td4([&]() {
        //打印直方图
        std::ofstream outFile;
        outFile.open(GetFilePath(args[0]) + "histogram.csv", std::ios::out); // 打开模式可省略
        outFile << "误差" << ',' << "数量" << std::endl;
        for (int i = 256; i < 512; ++i) {
            outFile << i - 512 << ',' << histogram[i] << std::endl;
        }
        for (int i = 0; i < 256; ++i) {
            outFile << i << ',' << histogram[i] << std::endl;
        }
        outFile.close();

        double meanSquareError = 0;
        double meanError = 0;
        int64_t weightsSum = 0;
        //计算平均值
        for (int i = 256; i < 512; ++i) {
            meanError += (double) (i - 512) * histogram[i];
            weightsSum += histogram[i];
        }
        for (int i = 0; i < 256; ++i) {
            meanError += (double) (i) * histogram[i];
            weightsSum += histogram[i];
        }
        meanError /= weightsSum;
        //计算方差
        for (int i = 256; i < 512; ++i) {
            meanSquareError += ((double) i - 512 - meanError) * ((double) i - 512 - meanError) * histogram[i];
        }
        for (int i = 0; i < 256; ++i) {
            meanSquareError += ((double) i - meanError) * ((double) i - meanError) * histogram[i];
        }
        meanSquareError /= weightsSum;
        meanSquareError = sqrt(meanSquareError);
        printf("meanSquareError = %f\r\n", meanSquareError);
        printf("meanError = %f\r\n", meanError);
    });
    td2.join();
    td3.join();
    td1.join();
    td4.join();
    std::system("pause");
    return 0;
}

int main(int argn, const char **argv) {
    if (argn == 1) {
        printf("文件名示例：\r\n");
        printf("**** 2^1us*****\r\n");
        printf("**** 2^1ms*****\r\n");
        printf("**** 1231us*****\r\n");
        printf("**** 123.4ms*****\r\n");
        printf("***中不能出现'^','us','ms',以及在关键词前留下空格\r\n");
        std::system("pause");
        return 0;
    }
    std::vector<std::string> args;
    if (strcmp(argv[1], "--test") == 0) {
        char *PIC_PATH = R"(C:\Users\17616\Desktop\本科毕设\实验数据记录\2021-5-18\G0_0-LD2-C0-I2-F4-0518\G0_0-LD2-C0-I2-F4-0518[%dus]png)";
        for (int i = 64; i < 8e6; i *= 2) {
            args.push_back(Format(PIC_PATH, i));
        }
    } else if (argn >= 4) {//保证有3张图以上
        for (int i = 1; i < argn; ++i) {
            args.emplace_back(argv[i]);
        }
    } else {
        return -1;
    }
    return m_main(args);
}


void calculate_roi(int rows, int cols) {
    float deltaR = (float) rows / BLOCK_SIZE;
    float deltaC = (float) cols / BLOCK_SIZE;
    float sumR = 0;
    float sumC = 0;
    std::vector<int> vr(BLOCK_SIZE + 1);
    std::vector<int> vc(BLOCK_SIZE + 1);
    vr[0] = 0;
    vc[0] = 0;
    vr[BLOCK_SIZE] = rows;
    vc[BLOCK_SIZE] = cols;
    for (int i = 1; i < BLOCK_SIZE; ++i) {
        sumR += deltaR;
        sumC += deltaC;
        vr[i] = sumR;
        vc[i] = sumC;
    }
    for (int i = 0; i < BLOCK_SIZE; ++i) {
        for (int j = 0; j < BLOCK_SIZE; ++j) {
            ROIs[i * BLOCK_SIZE + j] = cv::Rect(vc[j], vr[i],
                                                vc[j + 1] - vc[j], vr[i + 1] - vr[i]);
        }
    }
};
